/**
 * @file   Geometry.h
 * @author Robert Lie
 * @date   Fri Jun 23 15:44:30 2003
 * 
 * @brief  
 * 
 * 
 */

#ifndef _Geometry_h_
#define _Geometry_h_

#ifndef NULL
#define NULL 0
#endif // NULL

#include <stdarg.h>
#include <dlfcn.h>

#include <cmath>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <list>
#include <iterator>
#include <algorithm>

#include <base/exceptions.h>

#include "Miscellaneous.h"
#include "Migration.h"

AFEPACK_OPEN_NAMESPACE

template <int DIM> class Point;
class Geometry;
class GeometryBM;
template <int DIM, int DOW> class Mesh;
template <int DIM> class QuadratureInfo;
template <int DIM> class QuadratureInfoAdmin;
template <int DIM> class TemplateGeometry;
template <int DIM, int DOW> class SimplestMesh;

template <int DIM> Point<DIM> midpoint(const Point<DIM>&, const Point<DIM>&);
template <int DIM> double distance(const Point<DIM>&, const Point<DIM>&);
template <int DIM> Point<DIM> barycenter(const std::vector<Point<DIM> >&, const double * = NULL);
template <int DIM> Point<DIM> operator+(const Point<DIM>&, const Point<DIM>&);
template <int DIM> Point<DIM> operator-(const Point<DIM>&, const Point<DIM>&);
template <int DIM> std::istream& operator>>(std::istream&, Point<DIM>&);
template <int DIM> std::ostream& operator<<(std::ostream&, const Point<DIM>&);

std::istream& operator>>(std::istream&, Geometry&);
std::ostream& operator<<(std::ostream&, const Geometry&);
bool is_same(const Geometry&, const Geometry&);

template <int DIM> filtering_istream& operator>>(filtering_istream&, QuadratureInfo<DIM>&);
template <int DIM> std::ostream& operator<<(std::ostream&, const QuadratureInfo<DIM>&);

template <int DIM> filtering_istream& operator>>(filtering_istream&, QuadratureInfoAdmin<DIM>&);
template <int DIM> std::ostream& operator<<(std::ostream&, const QuadratureInfoAdmin<DIM>&);

typedef std::vector<Geometry>::iterator GeometryIterator;

template <int DIM, int DOW> std::istream& operator>>(std::istream&, Mesh<DIM,DOW>&);
template <int DIM, int DOW> std::ostream& operator<<(std::ostream&, const Mesh<DIM,DOW>&);

template <int DIM> filtering_istream& operator>>(filtering_istream&, TemplateGeometry<DIM>&);
template <int DIM> std::ostream& operator<<(std::ostream&, const TemplateGeometry<DIM>&);

///////////////////////////////////////////////////////////////////////////

/**
 * point of dimension \p DIM. The class provides facilities of operating
 * on the coordinate of a point.
 */
template <int DIM>
class Point
{
 public:
  enum { dim = DIM };
 private:
  double x[DIM]; /**< Coordinate of the point. */
 public:
  Point(); /**< Default constructor. */
  Point(const double *); /**< Constructor with data from a double array. */
  Point(const Point&); /**< Copy constructor. */
  Point(double, ...); /**< Constructor taking parameters as the entries of the coordinate. */
  ~Point(); /**< Destructor. */
 public:
  Point<DIM>& operator=(const Point<DIM>&); /**< Copy a point */
  operator const double *() const; /**< Casting to double pointer. */
  operator double *(); /**< Casting to double pointer. */
  const double& operator[](int) const; /**< Acess to entry. */
  double& operator[](int); /**< Acess to entry. */
  double length() const; /**< Length of the vector from the origin to the point. */
  Point<DIM>& operator+=(const Point<DIM>&);
  Point<DIM>& operator-=(const Point<DIM>&);
  Point<DIM>& operator*=(const double&);
  Point<DIM>& operator/=(const double&);
 public:
  friend Point<DIM> midpoint<>(const Point<DIM>&, const Point<DIM>&); /**< The middle point of two points. */
  friend double distance<>(const Point<DIM>&, const Point<DIM>&); /**< The middle point of two points. */
  friend Point<DIM> barycenter<>(const std::vector<Point<DIM> >&, const double *);
  friend Point<DIM> operator+ <>(const Point<DIM>&, const Point<DIM>&); /**< Add the coordinate of two point together. */
  friend Point<DIM> operator- <>(const Point<DIM>&, const Point<DIM>&); /**< Minus the coordinate. */
  friend std::istream& operator>> <>(std::istream&, Point<DIM>&); /**< Stream input. */
  friend std::ostream& operator<< <>(std::ostream&, const Point<DIM>&); /**< Stream output. */
};

namespace Migration{

  template <class BUF, int DIM>
    ostream<BUF>& operator<<(ostream<BUF>& os, const Point<DIM>& t) {
    os.encode_binary((void *)(&t[0]), sizeof(double)*DIM);
    return os;
  }
  template <class BUF, int DIM>
    istream<BUF>& operator>>(istream<BUF>& is, Point<DIM>& t) {
    is.decode_binary(&t[0], sizeof(double)*DIM);
    return is;
  }

} // namespace Migration

/**
 * the data to describe a geometry. We take the index of the vertices and
 * boundary geometries of geometry to describe the geometry. Information
 * else can be generated by iteration algorithm.
 */
class Geometry
{
 private:
  int ind; /**< Index of the geometry. */
  std::vector<int> vtx; /**< Index of vertices. */
  std::vector<int> bnd; /**< Index of boundary geometries. */
 public:
  Geometry(); /**< Default constructor. */
  Geometry(const Geometry&); /**< Copy constructor. */
  ~Geometry(); /**< Destructor. */
 public:
  Geometry& operator=(const Geometry&); /**< Copy a geometry. */
  int index() const;
  int& index();
  int n_vertex() const; /**< Number of vertices. */
  int n_boundary() const; /**< Number of boundary geometries. */
  const std::vector<int>& vertex() const; /**< The vertex index array. */
  std::vector<int>& vertex(); /**< The vertex index array. */
  const std::vector<int>& boundary() const; /**< The boundary geometry index array. */
  std::vector<int>& boundary(); /**< The boundary geometry index array. */
  int vertex(int) const; /**< An entry of the vertex index array. */
  int& vertex(int); /**< An entry of the vertex index array. */
  int boundary(int) const; /**< An entry of the boundary geometry index array. */
  int& boundary(int); /**< An entry of the boundary geometry index array. */
 public:
  friend bool isSame(const Geometry&, const Geometry&); /**< Judge if two geometries are the same. */
  template <int, int> friend class Mesh; /**< class \p Mesh can acess \p Geometry */
  friend std::istream& operator>>(std::istream&, Geometry&); /**< Stream input. */
  friend std::ostream& operator<<(std::ostream&, const Geometry&); /**< Stream output. */
};

/**
 * Geometry with boundary marker. This class provide an integer to differ the boundary
 * mark of different geometry. According the boundary marker of geometry, the package
 * will generate boundary marker for degrees of freedom.
 */
class GeometryBM : public Geometry
{
 public:
  typedef int bmark_t;
 private:
  bmark_t bm; /**< Boundary marker. */
 public:
  GeometryBM(); /**< Default constructor. */
  GeometryBM(const GeometryBM&); /**< Copy contructor. */
  ~GeometryBM(); /**< Destructor. */
 public:
  GeometryBM& operator=(const GeometryBM&); /**< Copy operator. */
  int boundaryMark() const; /**< Acess to the boundary marker. */
  int& boundaryMark(); /**< Acess to the boundary marker. */
 public:
  template <int, int> friend class Mesh;
  friend std::istream& operator>>(std::istream&, GeometryBM&); /**< Stream input. */
  friend std::ostream& operator<<(std::ostream&, const GeometryBM&); /**< Stream output. */
};

/**
 * The data structure of a mesh. The class \p Mesh administrate a set of points and
 * a set of geometries. The geometries are organized according its dimension and stored
 * in arrays. A lot of mechanism provided to retrieve information from the mesh.
 */
template <int DIM, int DOW=DIM>
  class Mesh
  {
  public:
  enum { dim = DIM, dow = DOW };
  typedef Point<DOW> point_t;
  public:
  typedef GeometryBM::bmark_t bmark_t;
  private:
  typedef Mesh<DIM,DOW> mesh_t;

  std::vector<point_t > pnt; /**< Point array of the mesh. */
  std::vector<std::vector<GeometryBM> > geo; /**< Geometries arrays of the mesh.
                                                The geometries in \p n dimension are in the \p n-th entry of the array, which
                                                is still an array. */
  public:
  Mesh(); /**< Default contructor. */
  Mesh(const mesh_t&); /**< Copy constructor. */
  virtual ~Mesh(); /**< Destructor. */
  public:
  mesh_t& operator=(const mesh_t&); /**< Copy operator. */
  unsigned int n_point() const; /**< Number of points in the mesh. */
  unsigned int n_geometry(int) const; /**< Number of geometries in certain dimension. */
  const std::vector<point_t >& point() const; /**< Point array. */
  std::vector<point_t >& point(); /**< Point array. */
  const point_t& point(int) const; /**< A certain point. */
  point_t& point(int); /**< A certain point. */
  const std::vector<std::vector<GeometryBM> >& geometry() const; /**< Geometries arrays. */
  std::vector<std::vector<GeometryBM> >& geometry(); /**< Geometries arrays. */
  const std::vector<GeometryBM>& geometry(int) const; /**< Geometries array in certain dimension. */
  std::vector<GeometryBM>& geometry(int); /**< Geometries array in certain dimension. */
  const GeometryBM& geometry(int, int) const; /**< Certain geometry in certain dimension. */
  GeometryBM& geometry(int, int); /**< Certain geometry in certain dimension. */
  bmark_t boundaryMark(int, int) const; /**< Boundary marker of certain geometry in certain dimension. */
  bmark_t& boundaryMark(int, int); /**< Boundary marker of certain geometry in certain dimension. */
  void renumerateElement(); /**< Renumerate the element of the mesh. This is a very simple
                               but efficient algorithm to decrease the band-width of the sparse system obtained. */
  void renumerateElementHSFC(void (*f)(const double*, double*)=NULL);
  void readData(const std::string&); /**< Read in data from a file in the internal data format. */
  void writeData(const std::string&) const; /**< Write data to a file in the internal data format. */
  void readData1d(const std::string&); /**< Read in 1 dimensional data from a file with only the node corrdinates. */
  virtual void writeEasyMesh(const std::string&) const {};
  virtual void writeTecplotData(const std::string&) const {};

  DeclException1(ExcMeshData, 
                 char *, 
                 << "Mesh data error: " << arg1);
  public:
  friend std::istream& operator>> <>(std::istream&, mesh_t&); /**< Stream input. */
  friend std::ostream& operator<< <>(std::ostream&, const mesh_t&); /**< Stream output. */
  friend filtering_istream& operator>> <>(filtering_istream&, TemplateGeometry<DIM>&);
  friend std::ostream& operator<< <>(std::ostream&, const TemplateGeometry<DIM>&);
  };

/**
 * Quadrature information of certain geometry. This class provide those data about the
 * coordinate and weight of certain numerical quadrature formula in different algebraic
 * accuracy. And methods to retrieve information from the data structure are provided. */
template <int DIM>
class QuadratureInfo
{
 public:
  enum { dim = DIM };
 private:
  int alg_acc; /**< Algebraic accuracy.*/
  std::vector<Point<DIM> > pnt; /**< The coordinate of quadrature point. */
  std::vector<double> wei; /**< Quadrature weight on the point. */
 public:
  QuadratureInfo(); /**< Default contructor. */
  QuadratureInfo(const QuadratureInfo<DIM>&); /**< Copy constructor. */
  ~QuadratureInfo(); /** Destructor. */
 public:
  QuadratureInfo<DIM>& operator=(const QuadratureInfo<DIM>&); /**< Copy operator. */
  int n_quadraturePoint() const; /**< Number of quadrature points. */
  int algebricAccuracy() const; /**< Algebraic accuracy. */
  int& algebricAccuracy(); /**< Algebraic accuracy. */
  const std::vector<Point<DIM> >& quadraturePoint() const; /**< Quadrature point array. */
  std::vector<Point<DIM> >& quadraturePoint(); /**< Quadrature point array. */
  const Point<DIM>& quadraturePoint(int) const; /**< Certain quadrature point. */
  Point<DIM>& quadraturePoint(int); /**< Certain quadrature point. */
  const std::vector<double>& weight() const; /**< Quadrature weight array. */
  std::vector<double>& weight(); /**< Quadrature weight array. */
  const double& weight(int) const; /**< Certain quadrature weight. */
  double& weight(int); /**< Certain quadrature weight. */
 public:
  friend filtering_istream& operator>> <>(filtering_istream&, QuadratureInfo<DIM>&); /**< Stream input. */
  friend std::ostream& operator<< <>(std::ostream&, const QuadratureInfo<DIM>&); /**< Stream output. */
};

/**
 * Administrator of a list of quadrature information. It provides retrieving facilities
 * to get quadrature information according required algebraic accuracy.
 */
template <int DIM>
class QuadratureInfoAdmin : public std::vector<QuadratureInfo<DIM> >
{
 public:
  enum { dim = DIM };
 private:
  std::vector<int> acc_tbl; /**< Algebraic accuracy table. */
 public:
  QuadratureInfoAdmin(); /**< Default contructor. */
  QuadratureInfoAdmin(const QuadratureInfoAdmin<DIM>&); /**< Copy Contructor. */
  ~QuadratureInfoAdmin(); /** Destructor. */
 public:
  QuadratureInfoAdmin<DIM>& operator=(const QuadratureInfoAdmin<DIM>&); /**< Copy operator. */
  const QuadratureInfo<DIM>& find(int) const; /**< Retrieve quadrature information according algebraic accuracy. */
  QuadratureInfo<DIM>& find(int); /**< Retrieve quadrature information according algebraic accuracy. */
 public:
  friend filtering_istream& operator>> <>(filtering_istream&, QuadratureInfoAdmin<DIM>&); /**< Stream input. */
  friend std::ostream& operator<< <>(std::ostream&, const QuadratureInfoAdmin<DIM>&); /**< Stream output. */
};

/**
 * Template geometry is the geometry information of a template element. A template
 * geometry is in fact a one-element mesh. A template geometry have the information
 * about how to calculate its volume. Such a function is stored in a shared library.
 * The user should provide such a shared library and tell this class about the file
 * name of the shared library and the function name to calculate the volume.
 */
template <int DIM>
class TemplateGeometry : public Mesh<DIM,DIM>
{
 public:
  enum { dim = DIM };
 private:
  std::string library_path;
  void * handle; /**< Pointer to object to open the shared library. Internel use only. */
  std::string library_name; /**< The file name of the shared library. */
  std::string volume_function_name; /**< The name of the function to calculate its volume. */
  double (*volume_function)(const double **); /**< The volume function pointer. */
  QuadratureInfoAdmin<DIM> quad_info; /**< The quadrature information on the geometry. */
 public:
  TemplateGeometry(); /**< Default contructor. */
  TemplateGeometry(const TemplateGeometry<DIM>&); /**< Copy constructor. */
  ~TemplateGeometry(); /** Destructor. */
 public:
  TemplateGeometry& operator=(const TemplateGeometry<DIM>&); /**< Copy operator. */
  void loadFunction(); /**< Load the function from the shared library. */
  void unloadFunction(); /**< Close the shared library. */
  const std::vector<Point<DIM> >& vertexArray() const; /**< Vertex array. */
  const QuadratureInfoAdmin<DIM>& quadratureInfo() const; /**< Quadrature information. */
  QuadratureInfoAdmin<DIM>& quadratureInfo(); /**< Quadrature information. */
  const QuadratureInfo<DIM>& findQuadratureInfo(int) const; /**< Retrieve quadrature information. */
  double volume() const; /**< Volume of the template geometry. */
  void readData(const std::string&); /**< Read in data from a file in certain file format. */
  void writeData(const std::string&) const; /**< Write out data to a file in certain file foramt. */

  DeclException1(ExcTemplateGeometryData, 
                 char *, 
                 << "Template geometry data error: " << arg1);
  DeclException1(ExcFileOpen, char *,
                 << "Can't open library "
                 << arg1);
  DeclException2(ExcLoadFunction, char *, char *,
                 << "Can't load function "
                 << arg1
                 << " from library "
                 << arg2);
	
 public:
  friend filtering_istream& operator>> <>(filtering_istream&, TemplateGeometry<DIM>&); /**< Stream input. */
  friend std::ostream& operator<< <>(std::ostream&, const TemplateGeometry<DIM>&); /**< Stream output. */
};

/**
 * Simplest mesh is a kind of mesh with only the points coordinate and the element
 * geometry information. Some grid generation program provide such kind of data format.
 * This class provides facilities to generate a mesh with internal data format from a
 * simplest mesh. This class will be helpful when the data provided by grid generation
 * program. Generally, a grid generation program will be sure to provide this class
 * required. Warning: to generate a mesh with internal data format will be really
 * time-consuming.
 */
template <int DIM, int DOW=DIM>
  class SimplestMesh
  {
  public:
  enum { dim = DIM, dow = DOW };
  protected:
  struct SimplestGeometry {
  public:
    int template_element; /**< Index of template element for the simplest geometry. */
    std::vector<int> vertex; /**< The vertex array of this element. */
  };
  std::vector<Point<DOW> > pnt; /**< Point array of the mesh. */
  std::vector<SimplestGeometry> ele; /**< Element array of the mesh. */
  std::vector<TemplateGeometry<DIM> > * tg; /**< Template element array for the mesh. */
  public:
  SimplestMesh() {}; /**< Default contructor. */
  virtual ~SimplestMesh() {}; /**< Destructor. */
  public:
  int n_point() const {return pnt.size();} /**< Number of points in the mesh. */
  int n_element() const {return ele.size();} /**< Number of elements in the mesh. */
  const std::vector<Point<DOW> >& point() const {return pnt;} /**< Point array. */
  std::vector<Point<DOW> >& point() {return pnt;} /**< Point array. */
  const Point<DOW>& point(int i) const {return pnt[i];} /**< Certain point. */
  Point<DOW>& point(int i) {return pnt[i];} /**< Certain point. */
  const std::vector<SimplestGeometry>& element() const {return ele;} /**< Element array. */
  std::vector<SimplestGeometry>& element() {return ele;} /**< Element array. */
  const SimplestGeometry& element(int i) const {return ele[i];} /**< Certain element. */
  SimplestGeometry& element(int i) {return ele[i];} /**< Certain element. */
  const std::vector<int>& elementVertex(int i) const {return ele[i].vertex;} /**< Vertex array of certain element. */
  std::vector<int>& elementVertex(int i) {return ele[i].vertex;} /**< Vertex array of certain element. */
  int elementVertex(int i, int j) const {return ele[i].vertex[j];} /**< Certain vertex of certain element. */
  int& elementVertex(int i, int j) {return ele[i].vertex[j];} /**< Certain vertex of certain element. */
  const std::vector<TemplateGeometry<DIM> >& templateGeometry() const {return *tg;} /**< Template element array. */
  std::vector<TemplateGeometry<DIM> >& templateGeometry() {return *tg;} /**< Template element array. */
  const TemplateGeometry<DIM>& templateGeometry(int i) const {return (*tg)[i];} /**< Certain template element. */
  TemplateGeometry<DIM>& templateGeometry(int i) {return (*tg)[i];} /**< Certain template element. */
  void setTemplateGeometry(std::vector<TemplateGeometry<DIM> >& t) {tg = &t;} /**< Set template array. */
  public:
  virtual void generateMesh(Mesh<DIM,DOW>&); /**< Generate mesh with internal data format. */
  };

AFEPACK_CLOSE_NAMESPACE

#endif //_Geometry_h_

/**
 * end of file
 * 
 */

